import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import seaborn as sns
import numpy as np
import plotly.express as px
from matplotlib import pyplot as plt
from scipy.stats.mstats import winsorize
from sklearn.decomposition import PCA
from sklearn.impute import KNNImputer
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.power import FTestAnovaPower
from interpret.glassbox import ExplainableBoostingRegressor
from interpret import show
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

coffee_df = pd.read_csv("CoffeeEconomies.csv")
"""Dropping meta-variables"""
coffee_df = coffee_df.drop(["i_cig", "i_xm", "i_xr", "i_outlier", "i_irr", "cor_exp", "statcap"], axis=1)

"""Lagging the dependent variable right away"""
coffee_df["coffee_val_1Y_lag"] = coffee_df.groupby('country')['coffee_val'].shift(-1)

"""Features-only dataframe"""
features_df = coffee_df.drop(["country", "coffee_val", "coffee_val_1Y_lag", "year"], axis=1)
feature_names = features_df.columns

"""Missing values imputation: for the purposes of VIF analysis and PCA!"""
features_imputed_df = pd.DataFrame(KNNImputer(n_neighbors=20, weights="distance").fit_transform(features_df), columns=feature_names)

"""Dropping redundant and highly correlated predictors"""
features_imputed_df = features_imputed_df.drop(["cgdpo", "cgdpe", "rgdpo", "rgdpe", "pop", "cn", "ck", 
                                                "rconna", "rdana", "cda", "ccon", "rgdpna", "ctfp", 
                                                "pl_con", "pl_da", "pl_gdpo", "rtfpna", "labsh", "rnna", 
                                                "cwtfp", "energy", "delta", "csh_c", "csh_g", "temp_chan", 
                                                "rwtfpna", "hc", "avh", "land_cov"], axis=1)

"""PCA for combating multicollinearity for price and political variables"""
features_imputed_df["price_index"] = PCA(n_components=1).fit_transform(features_imputed_df[["pl_g", "pl_m", "pl_x", "pl_n", "pl_i", "pl_k", "pl_c"]])
features_imputed_df.drop(columns=["pl_g", "pl_m", "pl_x", "pl_n", "pl_i", "pl_k", "pl_c"], inplace=True)

features_imputed_df["political_index"] = PCA(n_components=1).fit_transform(features_imputed_df[["rule_law", "govern_eff", "corrupt_control", "pol_stab", "voice", "reg_qual"]])
features_imputed_df.drop(columns=["rule_law", "govern_eff", "corrupt_control", "pol_stab", "voice", "reg_qual"], inplace=True)

"""VIF for assessing multicollinearity"""
vif = [variance_inflation_factor(features_imputed_df.values, i) for i in range(features_imputed_df.shape[1])]
vif_table = pd.DataFrame({"Predictor": features_imputed_df.columns, "VIF": vif}).sort_values('VIF', ascending=False)
print(vif_table) # everything's under 10

"""Correlation matrix: excluding highly correlated variables"""
correlations = features_imputed_df.corr()
sns.clustermap(correlations, xticklabels=features_imputed_df.columns, yticklabels=features_imputed_df.columns)
plt.show

          Predictor       VIF
5             csh_i  9.633083
2              rkna  8.926735
7             csh_m  8.761391
3               irr  8.506804
6             csh_x  5.663585
0          disaster  5.070893
1               emp  3.613716
9       price_index  2.372130
8             csh_r  1.505137
4                xr  1.310648
10  political_index  1.277693

<function matplotlib.pyplot.show(close=None, block=None)>

power = FTestAnovaPower().solve_power(effect_size=0.15, nobs=748, alpha=0.05, k_groups=11) # nobs=748 because of lag
print(power)  # 0.811 -> sample size is okay

0.8112649371872992

"""Creating the lagged version of the final dataframe + Winsorize 2%"""
df = pd.concat([coffee_df[["country", "year", "coffee_val_1Y_lag"]], features_imputed_df[["political_index", "price_index"]], 
                features_df[["xr", "rkna", "irr", "csh_m", "csh_r", "csh_i", "emp", "disaster", "csh_x"]]], axis=1) # concat indexes from the imputed df and other features from the non-imputed one
df = df.drop(df.index[-44:])
df['coffee_val_1Y_lag'] = winsorize(df['coffee_val_1Y_lag'], limits=(0, 0.01))
df["rkna/emp"] = df["rkna"] / df["emp"] # it will help to control for capital intensivity
df = df.drop(["rkna"], axis=1) # don't think we do need this anymore
df_reset = df.reset_index(drop=True)

print(df.columns)

"""Slope charts for Y"""
px.line(df, x='year', y="coffee_val_1Y_lag", color='country', title="Coffee Supply Chain Dynamics (2002-2018)").show()

"""Split around 2014, to prevent data leakage"""
train = df[df['year'] < 2014]
test = df[df['year'] >= 2014]

Y_train = train["coffee_val_1Y_lag"]
Y_test = test["coffee_val_1Y_lag"]
X_train = train.drop(["coffee_val_1Y_lag", "year"], axis=1)
X_test = test.drop(["coffee_val_1Y_lag", "year"], axis=1)
print(X_test.columns)

Index(['country', 'year', 'coffee_val_1Y_lag', 'political_index',
       'price_index', 'xr', 'irr', 'csh_m', 'csh_r', 'csh_i', 'emp',
       'disaster', 'csh_x', 'rkna/emp'],
      dtype='object')

Index(['country', 'political_index', 'price_index', 'xr', 'irr', 'csh_m',
       'csh_r', 'csh_i', 'emp', 'disaster', 'csh_x', 'rkna/emp'],
      dtype='object')

Модель:

Explainable Boosting Machine -- это часть семьи generalized additive models, которая позволяет количественно оценить важность признаков, а также автоматически выявить их взаимодействия. На вход идет массив признаков, на выходе имеем количественную оценку и ранжирование важностей признаков.

ebm = ExplainableBoostingRegressor(feature_types=["nominal", "continuous", "continuous", "continuous", "continuous", 
                                                  "continuous", "continuous", "continuous", "continuous", "continuous", "continuous", 
                                                  "continuous"])
ebm.fit(X_train, Y_train)

ExplainableBoostingRegressor(feature_types=['nominal', 'continuous',
                                            'continuous', 'continuous',
                                            'continuous', 'continuous',
                                            'continuous', 'continuous',
                                            'continuous', 'continuous',
                                            'continuous', 'continuous'])

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

"""Basic Metrics"""
Y_pred = ebm.predict(X_test)
mae = mean_absolute_error(Y_test, Y_pred)
rmse = np.sqrt(mean_squared_error(Y_test, Y_pred))
r2 = r2_score(Y_test, Y_pred)
print("MAE: ", mae, "RMSE: ", rmse, "R2: ", r2)

Y_train_pred = ebm.predict(X_train)
train_mae = mean_absolute_error(Y_train, Y_train_pred)
train_rmse = np.sqrt(mean_squared_error(Y_train, Y_train_pred))
train_r2 = r2_score(Y_train, Y_train_pred)
print("Train MAE: ", train_mae, "Train RMSE: ", train_rmse, "Train R2: ", train_r2)

"""Permutations"""
from sklearn.inspection import permutation_importance

importances = permutation_importance(ebm, X_test, Y_test, n_repeats=10, random_state=42)

# Sort features by importance
sorted_idx = importances.importances_mean.argsort()[::-1]

print("Permutation importances:")

for i in sorted_idx:
    print(f"{X_test.columns[i]:<20}: {importances.importances_mean[i]:.3f} ± {importances.importances_std[i]:.3f}")

"""Naive prediction baseline"""
from sklearn.dummy import DummyRegressor
dummy = DummyRegressor(strategy="mean").fit(X_train, Y_train)
dummy_mae = mean_absolute_error(Y_test, dummy.predict(X_test))
dummy_rmse = np.sqrt(mean_squared_error(Y_test, dummy.predict(X_test)))
print(f"Baseline MAE (mean prediction): {dummy_mae}") #  MAE is 2 times better than naive one; 
print(f"Baseline MAE (mean prediction): {dummy_rmse}")

"""Residuals plot"""
sns.histplot(Y_test - Y_pred, kde=True)
plt.axvline(0, color='r')
plt.show() # residuals show that some features necessary for good predictions post-2014 are lacking; need dummies for the trade wars maybe?

"""Actual vs. Prediction"""
plt.scatter(Y_test, Y_pred, alpha=0.5)
plt.plot([Y_test.min(), Y_test.max()], [Y_test.min(), Y_test.max()], 'r--')
plt.xlabel("Actual")
plt.ylabel("Predicted")
plt.show()

"""Residuals by country"""
results_df = pd.DataFrame({
    'country': test['country'],
    'residuals': Y_test-Y_pred,
    'predicted': Y_pred
})

sns.boxplot(x='country', y='residuals', data=results_df, orient="v", gap=0.5)
plt.xticks(rotation=90)
plt.axhline(0, color='r', linestyle='--')

MAE:  3601.5924990375306 RMSE:  6309.089139697104 R2:  0.9192037473927236
Train MAE:  791.9218802178859 Train RMSE:  1252.3189757112136 Train R2:  0.9929071001855435
Permutation importances:
emp                 : 1.761 ± 0.144
xr                  : 0.733 ± 0.054
rkna/emp            : 0.727 ± 0.077
country             : 0.093 ± 0.007
political_index     : 0.015 ± 0.004
irr                 : 0.006 ± 0.002
csh_m               : 0.004 ± 0.001
csh_x               : 0.002 ± 0.001
csh_i               : 0.001 ± 0.001
price_index         : 0.001 ± 0.001
disaster            : 0.001 ± 0.002
csh_r               : -0.002 ± 0.002
Baseline MAE (mean prediction): 9442.670041322312
Baseline MAE (mean prediction): 22447.20176428716

<matplotlib.lines.Line2D at 0x7881261bc050>

Метрики:

Были использованы метрики MAE, RMSE и R^2. Основными метриками для меня являлись MAE и R^2, поскольку MAE более устойчив перед выбросами, а R^2 показывает объясненную долю дисперсии, что важно, на мой взгляд при интерпретации результатов модели.

MAE и RMSE тестовой выборки были сравнены с наивным baseline предсказанием и показали хорошее приращение точности предсказаний (соотв. 3601 vs. 9442 и 6309 vs. 22447).

Все три метрики тестовой выборки были сравнены с метриками тренировочной выборки, чтобы оценить степень переобучения. На мой взгляд, "R^2 gap" в 7.3% позволяет судить о средней степени переобученности модели. Можно предпринять попытку бОльшей регуляризации модели. В целом, однако, текущий результат удовлетворителен (R^2 = 0.919 это очень хороший результат).

Были построены графики остатков, по которым можно сделать вывод, что остатки имеют выбросы в хвостах нормального распределения. Об этом позволяет судить и actual vs. prediction plot -- группа стран-олигополистов стабильно отклоняется по своим ошибкам от остальной массы поставщиков. Это следствие некоторой несбалансированности датасета, которая была ослаблена легкой винсоризацией. Ни лог-трансформация, ни square-root трансформация, ни фиктивные переменные для олигополистов -- ничего не помогло полностью устранить эту проблему, а местами результат оказался и намного хуже (лог-трансфорация сильно искажает датасет в данном случае). Был построен график ошибок по странам. Индонезия и Бразилия имеют наибольший разброс остатков, что подтверждает несбалансированность выборки. Тем не менее, модель имеет хорошие метрические показатели и из нее можно получить полезную информацию, поэтому текущая конфигурация проекта принимается как рабочая.

"""Let's cross-validate the model, but using TS-specific splits"""
from sklearn.model_selection import cross_val_score, TimeSeriesSplit
scores = cross_val_score(ebm, X_train, Y_train, cv=TimeSeriesSplit(), scoring='r2')
print(f"CV R²: {np.mean(scores):.3f} ± {np.std(scores):.3f}") # CV R²: 0.600 ± 0.238 -- acceptable for the dataset

CV R²: 0.600 ± 0.238

Кросс-валидация: Судить о модели по метрикам ее одной итерации недостаточно. Поэтому была проведена кросс-валидация, причем с учетом панельной структуры данных. Скоринг кросс-валидации осуществлялся по R^2 и результаты кросс-валидации оказались приемлимыми: CV R²: 0.600 ± 0.238.

Важность признаков: Проведена вторичная валидация результатов через permutation importances. Топовые фичи совпадают с финальными результатами EBM моделирования (они в ячейке ниже). Значит модель относительно стабильна. Какие это фичи (судя по финальной модели)? По убыванию:

1. Занятость
2. Интенсивность использования капитала
3. Валютный курс
4. Страна (фиктивная переменная)
5. Политический индекс
6. Процентная ставка
7. Ценовой индекс
8. Доля экспорта в ВВП
9. Доля импорта в ВВП
10. Доля остатков в торговом балансе в ВВП
11. Катаклизмы
12. Доля инвестиций в ВВП. Интерактивные фичи являются очень малозначимыми и ими можно пренебречь.

Исходя из топ-4 фичей можно сказать, что поставки кофе больше всего зависят от фундаментальных факторов производства -- труда и капитала на душу. Это соответствует неоклассической парадигме в экономике. Валютный курс и фиксированные страновые эффекты показывают значимость торговой политике в стоимостных цепочках -- не удивительно, учитывая сырьевой характер исследуемого блага.

Пользователи:

1. Государственные и негосударственные организации (при составлении индексов и рейтингов, и также при оценке рыночной ситуации очень полезно установить значимость факторов).

2. Трейдеры (как в виде фирм, так и частные) -- позволяет сфокусироваться при фундаментальном анализе на отдельных важнейших факторах.

3. Академики -- могут сформировать общее системное видение кофейной отрасли и связать отдельные важные признаки с современными теориями и парадигмами.

Ограничения:

1. Стабильные и последовательные данные, которые использовались в данной модели публикуются лишь раз в несколько лет (PwT). Данные ВБ меняются в методологии и могут вносить новые искажения в модель -- проблема.
2. Переменные кофейного рынка оказались очень коллинеарны, а данных немного. Поэтому пришлось ограничиться лишь небольшим наборов признаков. Не все теоретически интересные признаки поэтому получилось оценить.

# Let's interpret the result
show(ebm.explain_global())
show(ebm.explain_local(X_test, Y_test), 0)